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The dynamic response of a gas bubble entrapped in a cavity on the surface of a 
submerged solid subject to an acoustic field is investigated in the linear approxima- 
tion. We derive semi-analytical expressions for the resonance frequency, damping and 
interface shape of the bubble. For the liquid phase, we consider two limit cases: po- 
tential flow and unsteady Stokes flow. The oscillation frequency and interface shape 
are found to depend on two dimensionless parameters: the ratio of the gas stiffness 
to the surface tension stiffness, and the Ohnesorge number, representing the relative 
importance of viscous forces. We perform a parametric study and show, among oth- 
ers, that an increase in the gas pressure or a decrease in the surface tension leads to 
an increase in the resonance frequency until an asymptotic value is reached. 
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I. INTRODUCTION 



The volume pulsations of a gas pocket entrapped on a liquid-covered solid surface con- 
stitute a fundamental problem at the root of several applications in biology, microffuidics, 
cavitation and others. For example, the oscillatory flow induced by the pulsations causes 
a liquid motion which can be used to study the behavior of bacteria and cells under the 
action of shear 1-3 . In these conditions sonoporation of cell walls may occur, which would 
facilitate the uptake of drugs 1 and gene transfection . The flow induced by the oscillating 
gas pocket also induces mixing and streaming''. Under large-amplitude acoustic excitation, 
small gas bubbles issue from the gas entrapped in the cavities which greatly enhance sono- 
chemical reactions in a more controlled way than is possible in a conventional sonoreactor' . 
Microfabricated cavities on a silicon surface have been used to study controlled cavitation 
and bubble growth and collapse 8 ' 9 . 

Despite this wide range of applications, little is known about the dynamic response of a 
gas pocket on a submerged solid in an acoustic field. Miller and Neppiras et al. n recorded 
the acoustic response of multiple bubbles entrapped in a membrane. However, their size was 
not controlled, and no information about the response of the individual bubbles could be 
obtained. Rathgen et al. 12 studied the dynamics of periodic arrays of gas-filled micropores 
of controlled size on a solid surface. Using optical diffraction techniques, they were able 
to resolve in time, with a high accuracy, the nanometer-scale oscillations of the gas-liquid 
menisci driven by a sound field. However, they were unable to resolve the shape of the 
menisci in the course of the oscillations. 

Theoretical studies mainly focused on spherical bubbles in the bulk liquid whereas, for 
crevice bubbles, only approximate results exist. Miller & Nyborg 11 derived approximate 
expressions for the lowest resonance frequency and damping of a gas-filled pore on a solid 
surface under the assumption that the interface shape is parabolic. Their result is that the 
lowest resonance frequency /o of a cylindrical pore with radius a and depth h is approximately 
given by 



with k the polytropic index, A = a/h the aspect ratio of the pore, po the gas pressure 
when the interface is flat, a the surface tension coefficient, and p the liquid density. As an 



/o 




(1) 
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example, upon taking A = 1, this relation predicts natural frequencies of 176, 17.6, and 1.76 
kHz for gas pockets of air in water with equivalent spherical radii of 10, 100, and 1000 /jm. 
To obtain (1), an energy argument was used. Rathgen et al. 12 improved somewhat on this 
estimate by formulating the correct hydrodynamic problem for the liquid phase. However, 
they only solved the problem in an approximate way, retaining the parabolic approximation 
for the free-surface shape and also considering only the lowest resonance frequency. 

The purpose of the present work is to study the dynamics of the liquid-gas interface 
bounding the gas contained in a cavity at the surface of a solid in the linear approximation. 
We calculate the frequency, damping and surface shape of the linear normal modes of os- 
cillation of the system in the inviscid and viscous cases. In many situations the resonance 
frequency of the system is mainly determined by the inertia of the liquid, and hence can be 
calculated with sufficient accuracy from a potential flow model. We estimate the damping in 
two ways: from the potential flow solution by using a dissipation function method and, more 
accurately, by solving the time-dependent Stokes equations. The dynamics of the liquid-gas 
interface is found to depend on two dimensionless parameters: the ratio of the gas stiffness to 
the surface tension stiffness, and the Ohnesorge number, representing the viscous damping 
during one period of oscillation. 

II. PROBLEM FORMULATION 

Our aim is to describe the resonance frequency and interface shape of a gas bubble 
entrapped in a crevice. We model the crevice as a cavity with a circular mouth at the 
surface of an infinite solid submerged in an incompressible liquid (Fig. 1). The cavity has 
an aspect ratio 

A = ira 3 /V , (2) 

with a the mouth radius and Vq the cavity volume when the interface is flat. For a cylindrical 
cavity, A = a/h with h the depth of the cavity, but the results that follow hold for cavities of 
arbitrary shape. We introduce a cylindrical coordinate system (r, z), with the origin located 
on the axis of the cavity mouth at the level of the infinite solid plane. The liquid-gas interface 
is assumed to remain pinned at the circular edge of the cavity. The elevation of the free 
surface over the plane z = is described by T](r, t), and is assumed to be small compared to 
the radius of the cavity mouth, rj « a. 
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FIG. 1. A cavity with an entrapped gas bubble. The radius of the circular cavity mouth a is 
indicated, as well as the cylindrical coordinate system (r,z). The pressures in the gas and in the 
liquid are p g and p, respectively. The elevation of the free surface above the (z = 0)-plane is 
denoted by r](r, t). 

When the interface is perturbed, the compression of the gas and the surface tension of the 
liquid-air interface provide restoring forces and the interface will start to oscillate around its 
equilibrium position, which is assumed to be flat. The oscillating interface causes a velocity 
field u and pressure field p in the liquid phase (the gas flow is neglected here). To calculate 
the oscillation frequency, damping and shape of the interface we need to couple the normal 
stress in the liquid, derived from the velocity field, to the pressure in the gas, which results 
from the effect of surface tension and the gas compression and expansion. 

We use a linear theory, in which the time-dependence is assumed to be proportional to 
e^, so that r](r,t) = ^(r)e^, with complex eigenfrequency ( = iu — (3, where u = 2nf 
is the angular frequency and /3 the damping coefficient due to viscous dissipation in the 
liquid. This time-dependence is left implicit in the expressions that follow. The motion of 
the liquid-air interface, described by S(r, t) = z — r/(r, t) = 0, is coupled to the velocity field 
in the liquid via the kinematic condition 

- N + u-VS = 0. (3) 

Since we consider a small interface deformation rj, we neglect all terms which are of second 
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order and higher, which leads to a kinematic boundary condition in the form 

u z \ z=0 = (V, < r/a < 1. (4) 

The solid is impermeable for the liquid, and therefore we impose along the remainder of the 
(z = 0)-plane 

u z \ z=0 = 0, 1 < r/a < oo. (5) 
Furthermore, as there is no slip on the solid surface, 

u r \ z=0 = 0, 1 < r/a < oo. (6) 

On a clean liquid-gas interface, a no-shear-stress boundary condition applies, since the gas 
viscosity is much smaller than the liquid viscosity: 

rr Z =^( d ^ + d ^)=0, 0<r/a<l, (7) 



dz dr r 

with fi the dynamic viscosity. The presence of impurities modifies the interfacial behavior 
and may be modeled by a no-slip condition 15 

u r \ z=Q = 0, < r/a < 1. (8) 

Since the mixed boundary value problem (6), (7) leads to a mathematical problem which 
does not appear to be solvable, we impose the no-slip conditions (6), (8) on the entire surface 
in the following analysis. 

The coupling of the pressures in the liquid and the gas occurs via the dynamic boundary 
condition at z = 

Flu 

(9) 



p g = p + crC - 2 /i 



dz z=0 

with p g the pressure in the gas bubble and C = — (d^T] + d r r\/r) the curvature of the free 
surface; the last term on the right-hand side is the viscous normal stress on the interface. A 
general relation between the gas volume V and the gas pressure is given by the polytropic 
expression 

with k the polytropic index; k = 1 is applicable to isothermal conditions. When the interface 
is flat (77 = 0), V — Vq and p g = p , the ambient pressure in the liquid. By expanding (10) 
for small interface deformations, we find 

V 



Pg ~ PO 



1 - K 



Vo 



(11) 
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The gas volume can be found from the interface shape by integration: 

H~ 



V = V Q + 2vr / r]rdr = V ( 



1 + A- 



with 



H 



rr)(r)dr, 



(12) 



(13) 



proportional to the volume change of the gas due to the interface deformation. Combining 
(11) and (12), we obtain 

' (14) 



p g —po\l — k\ 
and the dynamic boundary condition (9) becomes 



Poll- k\- 



H 



p — a ( d 2 T r] + -d r i] J — 2 p 



du, 



dz 



(15) 



2=0 



In the following sections, the pressure in the liquid will be calculated in two limit cases: 
potential flow on one hand, as described in Section III, and unsteady Stokes flow on the other, 
in Section V. In addition, an estimate of the viscous damping is obtained using a modified 
potential flow model, where we calculate the dissipation in the bulk from the potential-flow 
solution (Section IV). Before we proceed, we introduce the following dimensionless quantities 

H 



V 



v = -• » = \/ — u, C = C, p 

which we will use from now on, thereby dropping the carets. 




a 3 p 



a 



-p, H 



(16) 



III. POTENTIAL FLOW 

To calculate the liquid pressure used in the dynamic boundary condition (15) we first 
neglect the influence of viscosity completely. Then, the flow is irrotational. As the resonance 
frequency of the system is mainly determined by the inertia in the liquid, it can be obtained 
from an inviscid flow model with a fair accuracy, as will be seen later. In the framework of 
potential flow theory, the no-slip condition (6) at the solid substrate in combination with 
either (7) or (8) cannot be enforced. 



A. Governing equations 



For irrotational flow, the velocity field can be written in terms of the potential as 
u = V0, which satisfies the Laplace equation V 2 = 0. At the liquid-gas interface <p has to 
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satisfy the kinematic boundary condition (4) 

u z \ z=Q =d z <P\ z=Q = Cv, 0<r<l, (17) 
and at the solid substrate the impermeability condition (5) 

u z \ z=0 = d z (j)\ z=Q = 0, 1 < r < oo. (18) 

We express the (dimensionless) pressure in the liquid using the linearized Bernoulli integral 

p= C<p, (19) 

a 

and, neglecting the viscous stress, find for the dynamic boundary condition at the gas-liquid 
interface (15) 

PH -d 2 rV -^d rV = (<j)\ z=0 , 0<r<l, (20) 



with 



P = — , (21) 

a 



the ratio of the gas stiffness to the surface tension stiffness. 

B. Reduction to an eigenvalue problem 

As shown in Appendix A, the Hankel transform can be used to reduce (17)-(20) to a set 
of dual integral equations in terms of <p. The solution of this system results in the following 
expression for the dynamic boundary condition (20) 

PH - dli] - -d r r] = -C 2 / v(k) J (kr)dk, (22) 
r Jo 

with 

v(k) = / 7]{s)sJ {ks)ds. (23) 
Jo 

Integration of (22) leads to the following eigenvalue problem to be solved for ( and rj 

V (r) = \pH{t 2 - 1) + C 2 V ~P- [Jo(k) - J (kr)} dk. (24) 

To find the solution of this integral equation we expand the interface deformation r\ into a 
Fourier-Bessel series 

oo 

V (r,t) = Y,Ck(t)J (j k r), (25) 



k=l 



7 



with jk denoting the feth zero of the Bessel function J . Substituting (25) into (24) and taking 
the inner product with rJ (j n r), we obtain the following generalized eigenvalue problem for 
the eigenfrequency ( (see Appendix A for details) 

2P V Jl0fc) Jl( - Jn) C k + l -C n f n Jl{j n ) = -C 2 f2ckJl(j k )Jl(jn)jkjnfUkJn), (26) 

Ik In 2 

k=l JK JTb k=l 

with f(jkJn) given by (A10). 



IV. WEAK VISCOUS EFFECTS 

In a real flow, viscous dissipation in both the bulk of the liquid and the boundary layer on 
the solid surface dampen the bubble oscillations. The ratio of dissipation in the boundary 



layer to dissipation in the bulk is given by a/ 5 , where 5 ~ y^/uj is the viscous boundary 
layer thickness. If we scale the angular frequency on the basis of the free bubble Minnaert 

1/2 

frequency, u ~ (^Po/p) /a, we find that the ratio of damping in the boundary layer to 
damping in the bulk is given by 

(27) 



upoa P 



pv 2 " AOh 2 ' 
with Oh the Ohnesorge number, defined by 



This dimensionless parameter is a measure of the damping during one period of oscillation. 
In case the bulk dissipation dominates, i.e. for smaller pits, we can use the potential flow 
solution to estimate the damping coefficient 16 . 

To describe the oscillations of the damped system, we use a Lagrangian formulation 
complemented by the Rayleigh dissipation function. We again express the interface shape 
in terms of the Fourier-Bessel series (25). The motion of the system is now given by 

^ f^i] _ _Oh— (29) 

dt \dc k J dc k dcu 

with Cfc = (ck, C = £k ~ £ P the Lagrangian, 1Z the Rayleigh dissipation function, defined as 

1Z = T>/2 with T> the rate of viscous dissipation in the liquid . To find an expression for 

C in terms of Ck, we calculate the kinetic energy and potential energy £ p of the system. 

The dimensionless kinetic energy of the liquid can be expressed as 18 

^ = \1<W = \\ A ^^ (3°) 



where n represents the unit surface normal directed out of the liquid. Due to the imper- 
meability condition (18), the integral (30) reduces to an integral over the bubble interface: 



£k = — 7r / 0(r, 0)d t r] rdr. (31) 
Jo 

The potential energy of the system can increase by an increase in area through the effect of 
surface tension or by a decrease in volume through compression of the gas 



£p= f dA- f ( Pg - Po )dV = vr / (d rV f rdr + \kPH\ 
Jaq Jvo Jo 1 

with P given by (21). The rate of viscous dissipation in potential flow reads 1 " 



(32) 



V = l h f (pL + **V dV = 20h / upnAA. (33) 
2 J v \dxk oxij J A dx k 



By integrating (33) over the surface we arrive at 



V = -871-Oh 



o 



dr 



d r T] rdr. (34) 



2=0 

Substituting (25) into (31), (32), and (34), we can express the kinetic energy, potential 
energy, and dissipation in terms of the degrees of freedom Ck and dk as 



oo oo 



£fc = 71 d k d njkjnJl(3k)Jl(jn)f(jk,jn)i (35) 

k I 

.00 co 00 

£p = ^£;Mm) +^pJ2J2fj Ji{jk)Ji{ji) ' (36) 

fe=l k=l 1=1 3k 3l 

00 00 

V = 87rOh^ ^ d kCnjkjnJl(jk)Jl(jn)9(jk,jn)- (37) 
k I 

with f(jk,jn) given by (A10), and g(jk,jn) given by (Bll). Substituting (35)-(37) into (29) 
and replacing dk by (ck, we obtain 

°° J ( ' ) J ( ' ) 1 

2P Y1 " ^~ 5_C * + ^ C nJl J l(in) = -^^cMj^JlUr^JkjnfUkJn) 

k=l 3k 3n 1 k=l 

K 

-4(Oh ^ C k Jl(jn)Jl(jk)Jnjkg(jn,jk)- (38) 
fc=l 

Note that putting Oh to zero in (38), i.e. neglecting the viscous dissipation, leads to exactly 
the same equation as derived before for potential flow, which reconfirms (26). 
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V. UNSTEADY STOKES FLOW 



In the previous sections potential flow was used to calculate the pressure and velocities 
in the liquid. On the solid substrate, however, the no-slip boundary condition (6) applies. 
Hence, a viscous boundary layer develops on the substrate, which cannot be accounted for in 
a potential flow model. Therefore, we repeat the calculation of the liquid pressure using an 
unsteady Stokes flow model. For small interface deformations 77 <C a, as is the case here, the 
nonlinear term of the Navier Stokes equations can be neglected with respect to the unsteady 
inertia term, and the unsteady Stokes equations describe the flow in the entire domain 11 '. 



A. Governing equations 

The unsteady Stokes equations in dimensionless form, with the dimensionless quantities 
as defined in (16), read 

(u = -Vp + OhVV (39) 

with Oh the Ohnesorge number defined in (28). We express the velocity in terms of a stream 
function \l/ as 

u = V x (-Veo) = — )e r + -— e z . 40 

\r J oz \r J r or 

Taking the curl of (39), we obtain 

£Q = OhV 2 ft, (41) 

with vorticity O = V x it = — V 2 fif/reg). As boundary conditions we have again the 
kinematic condition (4) and impermeability of the substrate (5). As mentioned before, we 
impose a no-slip condition on the solid (6) as well as on the bubble surface (8), to render 
the mathematical problem tractable. The dynamic boundary condition now reads 

PH -d?r)--d r r) = p\ z=0 , < r < 1, (42) 

with the pressure to be calculated from (39). Note that the normal viscous stress drops 
out from (42) as a consequence of (8) which, from the equation of continuity, implies that 
du z /dz = on z = 0. 
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B. Reduction to an eigenvalue problem 



The solution to (41) can again be expressed in terms of the Hankel transform. Then, the 
eigenvalue problem to be solved for (, r\ is given by (details of the calculation can be found 
in Appendix B) 



77(f) = ^PH(r 2 -l)- 



v{k) 



[Jo(fc) - Mkr)] \( 2 + COh (k 2 + k^k 2 + C/Oh 



dfc, < r < 1. 
(43) 



Again, we expand the interface deformation 77 into the Fourier-Bessel series (25), and obtain 
(see Appendix B for details) 



oo 1 
k=1 JkJn * 



CkjkjnJl(jk)Jl(jn) 



k=l 



J 2 o(s) 



(f k -S 2 )(3 2 n-S 2 ) 



C 2 + COh \s 2 + sy/s 2 + C/Oh] ) ds. (44) 



VI. NUMERICAL SOLUTION METHOD 



To find the resonance frequency and interface shape in the potential flow model, the 
generalized eigenvalue problem (26) is truncated to K terms and solved numerically with 
Mathematica 8 (Wolfram Research) for ( and see Appendix A. We studied the con- 
vergence of the sum (25) for the first three modes by taking up to K = 100 terms into 
account for P = to P = 500. We found that the system converges rapidly: for the lowest 
mode K = 3 was already sufficient for accurate reconstruction of the interface shape. For 
higher modes, the matrix size increases because more Bessel functions are required to de- 
scribe the interface shape: for mode 3, we used K = 6. The larger the matrix, the more 
eigenfrequencies can be calculated. 

The generalized eigenvalue problem (44) for the Stokes flow model has to be solved 
iteratively, due to the complexity of the integral. To this end, we split the integral into three 
parts, so that the equation to be solved becomes 

oo - 

2p J2 Jl(jfc ) Jl(jn ) + - Cn j 2 n J 2 (j n ) = 

Ikln 2 
k=K JKJn 

oo 

C kjkjnJl(jk)Jl(jn) [CfUkJn) + (Oh {g(jkjn) + h(j k , j n , (, Oh)}] , (45) 

fe=l 
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with f(jk,jn) given by (A10), g{jk,jn) by (Bll), and h(j k ,j n ) by (B12). To reduce (B13) 
to a generalized eigenvalue problem that can be solved with Mathematica, we write d = £c, 
so that the resulting system becomes 

' ^ / C ' +1) ^ C (i+1) ( ' II " I , (46) 





2PB+|C Oh[G + H(C (i) ,Oh 

with % the iteration number, c a vector with elements cj., the zero matrix , I the unit matrix, 
A given by (A12), B by (A13), C by (A14), G by (B14), and H by (B15). The matrices A, 
G, and H correspond to the three terms in which the integral in (44) is decomposed. The 
generalized eigenvalue problem (46) is solved iteratively with the modified potential-flow 
solution used as initial guess in matrix H. With this initial guess, we evaluate the 
integral (B12) numerically, which we then use to solve (46). This step permits an improved 
estimate of the eigenvalue, which is substituted again into the matrix H. This procedure 
is repeated until convergence is reached, i.e. until the difference in both frequency and 
damping between the current and the previous iteration is less than 0.01% of the current 
result. Again, we investigated the influence of the matrix size on the calculation of the 
eigenfrequency and interface shape for the lowest three modes by taking a system size up 
to K = 16. For the lowest two modes, K = 3 was sufficiently accurate to calculate the 
eigenfrequency, whereas for mode 3, K = 6 was used. Convergence of (46) was achieved 
after 4 iterations for the three lowest modes; see Appendix B for details. To further check 
the convergence of the solutions obtained, we used the potential flow solution as initial guess, 
and slowly increased the Ohnesorge number from to Oh s = 0.0303. Using this method, 
we obtained the same results for the resonance frequency and the damping as by starting 
directly at Oh s with the modified potential flow solution as initial guess. 



VII. RESULTS 

In the generalized eigenvalue problem for (, rj, only two dimensionless parameters appear: 
the ratio of the gas stiffness to the surface tension stiffness, P (21), and the Ohnesorge 
number Oh (28). For a gas pocket with A = 1 in a 15-//m cylindrical micropit submerged in 
water under standard conditions, the corresponding values of the dimensionless groups are 
P = P s = 20.5 and Oh=Oh s =0.0303. In this example, we find the dimensional resonance 
frequency for the first three modes in potential flow to be 121, 274, and 556 kHz, respectively. 
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FIG. 2. (Color online) The interface shape as calculated from the potential flow model for (a) 
P = 20.5, mode (black), mode 1 (blue online), and mode 2 (red online), (b) The interface shape 
for mode with P = (black), 40 (blue online), 80 (red online), 160 (orange online), 320 (green 
online). Note that for P > 80, mode has an extra node in addition to the one at the rim of the 
pit. 



The first result is not very different from the frequency estimate (1) by Miller & Nyborg 11 
which is 151 kHz for the lowest mode of a 15-/xm pit. The interface shape for the first three 
modes is depicted in Figure 2a. One can see that, with these parameter values, the largest 
contribution to the interface shape of mode comes from the first term in the Fourier-Bessel 
series (25), whereas for mode 1 the largest contribution comes from the second term, etc. 

Figures 3, 4, and 5 show how the resonance frequency and damping coefficient depend on 
the two dimensionless parameters P and Oh. As expected, both the damping and frequency 
increase with the mode number with the result that, after a generic initial perturbation, 
the bubble will oscillate the longest at its fundamental resonance frequency whereas higher 
frequencies dampen out earlier. The difference in resonance frequency between the potential 
flow (PF), modified potential flow (mPF), and Stokes flow (SF) models is very small which 
means that, in the parameter range of interest, the resonance frequency is mainly determined 
by inertia and can be obtained from the potential flow model with sufficient accuracy. 

In Fig. 3a, the graph of the frequency / = u/2ir versus P for mode shows that the 
resonance frequency first increases with P, until it levels off to a dimensionless value / = 1.72 
for P > 200, approximately. In the approximate solution (1) by Miller 11 such a plateau is 
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FIG. 3. (a) The dimensionless frequency / = oj/2tt of mode versus the dimensionless parameter 
P, defined in (21), for potential flow (black, solid), modified potential flow (black, dotted) and 
Stokes flow (black, dashed) with Oh=0.0303. The results for potential flow and modified potential 
flow nearly overlap. For comparison, the result of Miller & Nyborg (1) is also shown (dash- 
dotted), (b) Dimensionless damping j3 versus P for the modified potential flow and Stokes flow 
models. 



not observed. Figure 4a shows graphs of / versus P for modes 1 and 2. Here, a similar 
increase in / with P is observed, but the plateau is reached at larger values of P. Initially, 
the resonance frequency increases with P because at larger P it becomes more difficult to 
change the volume of the gas, and hence the system becomes stiffer, which leads to a higher 
resonance frequency. The reason for the occurrence of a plateau in the frequency lies in the 
increasing stiffness of the gas. As P increases and the system becomes stiffer, it becomes 
more difficult to decrease the gas volume change H, as defined in (13), and the system 
responds by increasing the area of the interface instead (see Fig. 2a). To further decrease 
the volume change, at some point an extra node has to appear in the interface shape, as 
can be seen in Fig. 2a for P > 80. This node is pushed towards the axis of the pit as 
P is increased further. In this way, the interface area increases more and more, and the 
net volume change due to the surface elevation eventually tends to zero, which means that 
V — > Vq. Figure 6 shows that this decrease in the amplitude of the volume oscillations 
occurs in such a way that the product PH tends to a constant value. Hence, the resonance 
frequency levels off, and the system eventually oscillates with a fixed interface shape. The 
interface shapes corresponding to modes 1 and 2 are depicted in Fig. 7. For mode 1, the 
extra node appears around P ~ 500, whereas for mode 2, it will occur at a larger P. 
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FIG. 4. (Color online) (a) The dimensionless frequency / of mode 1 (blue online), and mode 2 (red 
online), versus the dimensionless parameter P, defined in (21), for potential flow (solid), modified 
potential flow (dotted) and Stokes flow (dashed) with Oh=0.0303. A plateau in the frequencies 
similar to that for mode is observed, but it occurs at larger P (not shown in the figure for mode 
2). (b) Dimensionless damping /3 versus P for the different cases. 
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FIG. 5. (Color online) (a) The dimensionless frequency / of mode versus the dimensionless 
parameter P, defined in (21), for modified potential flow (dotted), and Stokes flow (dashed) with 
Oh=0. 01515 (blue online), Oh=0.0303 (black), and Oh=0.0606 (red online), (b) Dimensionless 
damping /3 versus P for the different cases. 



The parameter P also has an effect on the damping coefficient as shown in Figs. 3b and 
4b. For mode the damping coefficient increases with P, until a final plateau is reached. For 
mode 1, however, a minimum is observed and the plateau is reached for larger P. For mode 
2, the damping decreases and reaches a minimum beyond the maximum value of P shown 
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FIG. 6. (Color online) The relative stiffness of the gas pocket times its volume change, PH, versus 
P, for mode (solid) and mode 1 (dashed). Initially PH increases, until a maximum is reached, 
then it tends to a finite value. 
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FIG. 7. (Color online) The interface shape as calculated from the potential flow model with 
Oh=0.0303 for P = (black), P = 250 (blue online), and P = 500 (red online) for (a) mode 1 and 
(b) mode 2. Note that, for P = 500, mode 1 has developed an extra node. 



in the graph. The presence of a minimum could be explained as follows: as the stiffness 
of the gas increases with P, the relative volume change decreases, which leads to a smaller 
liquid displacement and viscous energy dissipation. However, as the frequency increases, the 
damping increases as well. These two effects compete, and give rise to a minimum in the 
damping coefficient. The plateau is reached at larger P, when the product PH tends to a 
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finite value (see also Fig. 6). The difference in damping coefficients between the mPF and 
SF models is larger for the higher modes. One also observes that, for some values of P, the 
mPF damping is even larger than the SF damping. This behavior is due to the difference in 
velocity profiles between the mPF and SF, and is also known to occur for shape oscillations 
of drops and bubbles 19 . 

Whereas P has a large influence on the resonance frequency of the pit, the influence of 
Oh is only very small as shown in Fig. 5. The influence of the Ohnesorge number on the 
damping coefficient is of course large. 



VIII. CONCLUSION 

The resonance frequency, damping and interface shape of a gas pocket entrapped on the 
surface of an submerged solid have been calculated. To describe the hydrodynamic prob- 
lem in the liquid domain, both a potential and an unsteady Stokes flow model have been 
used. The potential flow model gives a reliable prediction of the resonance frequency of 
the gas pocket, which is mainly determined by inertia in the liquid. To derive an estimate 
for the damping of the oscillations, the bulk dissipation was calculated from the potential 
flow model. A more accurate prediction of the damping was derived based on the unsteady 
Stokes flow model, which is valid throughout entire domain and therefore includes the con- 
tributions of both the boundary layer and the bulk. However, the Stokes flow results will 
overestimate the real damping somewhat in the case that the liquid-gas interface is clean: 
in the method described here, a no-slip condition on the free surface was used, which leads 
to some additional dissipation. 

The resonance frequency, damping and interface shape of an entrapped gas pocket depend 
on two dimensionless numbers: the ratio P of the gas stiffness to the surface tension stiffness, 
defined in (21), and the Ohnesorge number (28), which represents the relative importance 
of viscous forces. In general, the resonance frequency increases with increasing gas stiffness. 
However, an unexpected feature of our results is that, when the volume stiffness of the gas 
pocket greatly exceeds the surface stiffness, the normal modes develop an extra node, and 
the resonance frequency tends to an asymptotic value. 
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Appendix A: Potential flow solution 

To obtain a solution for the velocity potential we set 

POO 

d z( f)= / m(k)J (kr)e~ kz dk (Al) 
Jo 

and, upon integration, find that 

poo 

§(k)J (kr)e- kz dk. (A2) 



'o J o 
and, therefore 



Using the boundary conditions (17) and (18) and the orthogonality relation for Bessel func- 
tions we obtain 

poo pi 

kr$(k)J (kr)J (hr)drdk = / $(k)8(h - k)dk = ( / r](r)rJ (hr)dr, (A3) 

Jo Jo 

$(h) = ( [ v(r)rJ (hr)dr = (v{h). (A4) 
Jo 

Thus 

POO 

0(r,O) = -C/ v(k)J (kr)dk. (A5) 
Jo 

Substituting this result into (20) we find (22). 

The next step is to express the interface deformation 77 in terms of the Fourier-Bessel 
series (25), to obtain 

v(s) = f> fr J (sr)J (j k r)dr = V c k 3k J ^ J ^ , (A6) 

00 „x 00 

H = 2Y,c k rJ (j k r)dr = 2Y,-JiUk)- (A7) 
k=i J° k=i ^ k 

Substitution into (24) leads to 

E ckJoUkr) = \p f ^ )( ^ -i)+c 2 E c kjMjk ) r [Jo % Jo %l Ms) ds. 

k=l 1 k=l Jk k=l J ° [j k - 3 ) S 

(A8) 
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To obtain an equation for each of the unknowns, we multiply (A8) by rJo(j n r) and integrate 
between and 1, to find 

2P^—^-J l (i k )J l {] n ) + -CnilJKin) = -C 2 ^2ckjkJiUk) I . 2 °_ 2 ds / r J (sr) J (; n r) dr 

fc=1 JkJn * k=1 JO Jk S JO 



'0 Jfc S JO 

-C^CkJl (jk)jkjn J 1 (jn ) / ( jfc , jn ) , ( A9 ) 



fe=l 

with 

JiFil W 1 ' 1; I' 1' 1' ~ 1; I' I' I; > if ^ n ' (A10) 

^ 2 F 3 (2,2;|,|,|;-ji), if * = n, 

where 2-P3 is the hypergeometric function ". After truncation of the Fourier-Bessel series 
(25) to K terms, we can express the eigenvalue problem in matrix form as 

C 2 A + 2PB + ^ c = 0, (All) 

with A a K x A^-matrix with coefficients 

A kn = Jl(jk)Jl(jn)jkjnf(jk,jn), (A12) 

B a A" x A-matrix with coefficients 

R Jl(jk)Jl(jn) , A1 -v 

JkJn 

and C a K x K diagonal matrix with coefficients 

C kn = jlJi(jk)hn, (A14) 

and c a A-array with coefficients c\~. The generalized eigenvalue problem (All) is the solved 
with Mathematica 8 (Wolfram Research). 

Appendix B: Stokes flow solution 

The general solution to (41) in terms of the Hankel transform reads 

n(r,z)= [ G(k)Jx(kr)e- z ^ k2+c/0h dk, (Bl) 
Jo 
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with G to be determined from the boundary conditions (4), (5), (6), and (8). From (Bl) we 
can determine an expression for stream function \I/ using 

n = -V 2 (X) . (B2) 

The homogeneous solution of (B2) reads 

POO 

* fc ( r , z) = r / F{k)J l (kr)e- kz dk, (B3) 
Jo 

with F to be determined from the boundary conditions. The particular solution can be 
found from (41): 

n = -V* (£e.) = ^ Vfl, (B4) 

and hence 

Oh f 00 / 

%{r,z) = -r^- I G(k)Ji(kr)e- z V k2 +</° h dk. (B5) 
C Jo 

Once we know $ = ^ + we can find expressions for the velocity field 

/•OO Ql POO 

/ kF{k)Ji{kr)e- kz dk-— ^W^CjO^G{k)Ji{kr)e' z y k2+cl ° h dk, 
Jo C Jo 



u r {r,z) 
uJr, z) 



POO Ql POO 

/ kF(k)J (kr)e- kz dk- — / fcG(fc) J (fcr)e- z V fe2 +f/° h dfc. (B6) 



Expressions for F and G can now be obtained from (4), (5), (6), and 

F(k) - ^G(k) = (v(k), 
Oh 



kF(k) - — G(A)\/FTc70h = 0, 
with v(k) given by (23), which give 

G(k) = v(k) (k 2 + k^/k 2 + C/Oh) , (B7) 



F(ife) =v(A;) 



C + ^^ + ^v/p + c/Oh 



(B8) 



Oh 

The pressure in the liquid can now be obtained from (39). Using (40) and (41) one finds 
that 

Vp = -C (V x ^egj , (B9) 
and hence the liquid pressure reads 

POO 

p(r,z) = ( F(k)J (kr)e- kz dk. (BIO) 
Jo 
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This results evaluated at z = is then substituted into (42) to find (43) by integration. 
Substitution of the Fourier-Bessel series (25) results in (44). 

Due to the complexity of the integral in (44), the system has to be solved iteratively. To 
this end, we split the integral into parts. The resulting equation is given by (45) with / 
given by (A10), g by 

f^K 2 ^(i,i;i4^ Mn '(pii) 

-i 2 F 3 (l,2;|,|,|;-j|), if k = n. 

and h by 

h(j k ,j n ,C,Oh)=J o u , - _ g2) y/s* + C/Ohds. (B12) 

After truncation of the Fourier-Bessel series to K terms, the resulting equation (45) in 
matrix form becomes 

| (C (m) ) 2 A + 2PB + \ C + C (m) Oh [G + H(C (i) , Oh)] | c< i+1 ) = 0, (B13) 

with i the iteration number, A given by (A12), B given by (A13), and C given by (A14); G 
is a K x if-matrix with coefficients 

G kn = Jl(jk)Jl(jn)jkj n g(jk,jn), (B14) 

and H is a K x i^-matrix with coefficients 

H kn = Jl(jk)Jl(jn)jkjnh(j k ,jn, C W , Oh). (B15) 

Figure 8 shows that for mode 0, convergence is reached within 4 iterations, irrespective of 
the matrix size. 
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